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*vj This problem arises, for instance, in some homeland security applications. The goal 

C^^ is to reach the signal-to- noise ratio (SNR) levels on the order of 10^^. A Bayesian 

y~~^ approach to this problem is implemented in 2D. The method allows inference not 
T — I 

L| only about the existence of the source, but also about its location. We derive Bayes 

k>< factors for model selection and estimation of location based on Markov Chain Monte 

^ Carlo (MCMC) simulation. A simulation study shows that with sufficiently high total 

emission level, our method can effectively locate the source. 
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1. Introduction 

We consider the problem of detecting existence of a low emission radiating source inside a 
volume, in the presence of a strong random background. One can easily imagine possible 
applications of such detection, for instance to homeland security. We are interested 
in the situation when about 99.9% of the total detections come from the background 
particle^P] and from the particles emitted by the source that have been scattered (we 
will consider the latter as a part of the background). In other words, only about 1% of 
detected hits are by the ballistic particles coming from the source. Although medical 
emission tomographic imaging faces similar problems (e.g. p^), the overwhelming level of 
noise that has just been mentioned would be considered impossible to handle there. So, 
how can one attack this problem? Although there is probably no general solution, in the 
applications we have in mind, the radiating source, if present, would be significantly (on 
the order of hundred times) geometrically smaller than the whole object. As is explained 
in [I] , this, and the availability of detectors determining direction of an incoming particle, 
make detection conceivable under appropriate conditions. In this text, we consider the 
2D problem. Unlike 11], where more analytic techniques are considered, we propose a 
Bayesian method, which allows inference about the existence and location of the possible 
source. 

The problem can be stated as follows. One is interested in certain type of particles 
(say, 7-photons or neutrons, although the type of particles is irrelevant for our purpose). 
Suppose that the observed area belongs to the unit disk D in the /i-/2-plane (see Figure 
[I|. Detectors, placed around the object, are assumed to be able to determine the 
linear trajectory of each incoming particle. It is assumed that detectors surround the 
object in such a way that any escaping particle hits a detector (this assumption could 
in principle be weakened). Most of (or maybe all) detected particles are coming from 

correspondence author will be X.Xun. 

I? Most of the background particles are not emitted inside the object, but rather are present in the 

surrounding environment (e.g., cosmic rays). 
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Figure 1. A direction sensitive (e.g., collimated or Compton type) detector 
determines tlie normal parameters (9, S) of the trajectory of the incoming particle. 
The detected particles might be either emitted from the source or coming from random 

background. 
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a random background (and, in particular, are not emitted inside the object). Besides 
the background emission, a small (in comparison with the total object's size) and low 
emission source might be present. Many of the particles emitted by the source will 
be scattered, and only a small number of them will reach the detectors unscattered 
(ballistic). The goal is to detect the presence of such an object, if the emission is 
dominated by the background, e.g., such that the ballistic particles coming from a 
possible source could account to about 0.1% of the total emission. 

The set-up is illustrated in Figure [TJ The trajectory of a particle that hits the 
detector can be identified by its normal coordinates {9, S), and thus we assume that the 
detector provides the values (6*, S) for each hit. We expect that a radiation source 
of a small radius d might be present, in which case we denote its location point 
as L = {ii,i2)- If a particle is emitted from this source and reaches the detectors 
ballistically (unscattered), then 6 and S satisfy the inequality 



\ii cos(^) + £2 sin(0) -S\<d. 



(1) 



Most particles from the random background normally will not satisfy this condition, but 
a small portion might. The idea is that if a source is present, ballistic particles coming 
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from it might lead to a statistically significant increase in the number of trajectories 
satisfying ([T]), and thus to detection. Under appropriate conditions (geometrically 
sufficiently small source and sufficiently large total count in the sample), this happens 
to be the case (see the discussion in (T|). 

An outline of the paper is as follows. Section |2] introduces the candidate models 
describing the situations without a source and with a source, respectively. Section |3] 
explains the calculation of Bayes factors for determining the presence of a source, as 
well as the computational details of our Markov Chain Monte Carlo (MCMC) algorithm. 
Readers interested mostly in the viability of the approach can skip the subsections, 
going directly to Section |4] with the results of a simulation study, where the method is 
examined for various levels of source emission rate. These results confirm the possibility 
of detection. The article ends with conclusions and remarks. 

2. Models 

Suppose that the direction sensitive detectors registered hits by n particles and recorded 
the corresponding normal coordinates {6i,Si) for i = l,---,n of their incoming 
directions. We denote by 6i the (unobserved) indicator that the i^^ particle is coming 
ballistically from the suspected source. In other words, 5j = 1 if the z*'* particle comes 
from the source. Otherwise, 6i = 0. 

If there is no source present, then pr((5j = 1) = 0. If there is a source, we assume 
that the 5j's are independently and identically distributed according to the Bernoulli(p) 
law. This covers also the possible absence of a source, in which case p = 0. 

Our plan is to develop statistical models for each of these cases, and then decide, 
based on the collected data and the value of the corresponding Bayes factor, which 
model fits better the collected data. 

We would like to point out that this is not expected to be a simple problem, since 
it involves inference as to whether a non-negative parameter p takes its boundary value 
p = 0. Even in simple variance components models, frequentist boundary value testing 
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is a difficult matter, see for example [sl. 

2.1. The Model without a Source 

When there is no source (we call this model Mi), all hits at the detectors come from 
the random background and thus Si = for alH = 1, . . . , n. We will assume in this text 
that the random background is isotropic and uniform. In other words, the angle 6 and 
the distance S from the origin of the trajectory are uniformly distributed: 

[9i\5i = 0] =Uniform(0,27r); 

(2) 
[5.|5^ = 0] = Uniform(-l,l). 

Notice that particles having trajectories with 15*1 > 1 do not get detected and thus do 

not enter the model. 

2.2. The Model with a Source 

If a source exists (model M2), then pr(5j = 1) = p > 0. If the particle comes from 
the background, then Si = and relations ^ still hold. Assuming that the source in 
question is isotropic and uniform, when (5j = 1, we have the following distributions: 

[9i\6i = 1] =Uniform(0,27r); 

(3) 
[SilOi, (£1,4), Si = l]= ii cos(^i) + £2 sin(^i) + Uniform(-rf, d). 

Our goal thus is to choose between the models Mi and M2, based on the measured data. 

We show in the following subsection the priors, likelihood and the posterior distribution 

associated with model M2. 

2.2.1. Likelihood and Posterior of Model with a Source For the model M2 with a source, 
see (|3]), parameters of interest are = (p, £i,£2)- Priors are assigned as follows: 

L = (£1, £2) = Uniform{£2 + £2 < i|^ 

i.e., a priori the source can be located anywhere with equal probability. Also, we allow 
p to have the uniform discrete distribution on the set {0 < pi, ■ ■ ■ ,Ph}, which contains 
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h values equally spaced on the interval [pi,p/i]. This set is chosen based on a priori 
estimates of the possible emission strength of the source. 



For the implementation of the MCMC algorithm, it is convenient (see Section 3.1.2) 
to use the polar coordinates {r,u) of the center L: 

£i = rcosu, £2 = rsmu. 
Denoting the new parameterization hy ip = {p,r,u), the prior /(^|M2) is 

/(^iMs) = f{p,r,u\M2) = h-^n-^rl{0 <r < 1)1(0 <u< 27r), 
where !(■) is the indicator function. 



Let Yi = {6i, Si) denote the i observation (i.e., the normal coordinates of the 



■th 



trajectory (at the detector) of the i detected particle) and, as before, let 5, be the 
(unobserved) indicator associated with it. It will be convenient to introduce vectors 
Y = (Fi, . . . , Yn) and 6 = {6i, . . . ,6n)- Then the likelihood function is 

fiY\p,r,u,M2) (4) 

n 

= n(47r)"^ [pd~H {\r cos{u) cosiOi) + r sin(M) sin(6'i) - Si\ < d} + 1 - p] 

i=l 

n 

= (47r)""n [pd~^l{\rcos{u-9i)-Si\ < d} + 1 - p] 

i=l 

= (47r)-"(l - pY-\pd-^ + 1 - pY, 

where J = X]r=iHkcos('U — 9i) — Si\ < d} counts the total number of particles whose 
incoming trajectories at the detectors pass near the location {r,u). 
Given the above prior and likelihood, the posterior is 

f{^\Y, M2) = fip, r, u\Y, M2) (X r(l - p)"-^(prf-i + 1 - p)' , (5) 

where p G {pi, ■ ■ ■ ,Ph}, r G [0, 1] and u G [0, 27r]. 

3. Model Selection via Bayes Factors 

Let pr(Mj) be the prior probability of model Mj and pt{Y\Mj) be the marginal 
distribution of the data, given model Mj, where j = 1, 2. We also denote by pr(Mj|F) 
the posterior probability of the model Mj. 
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The parameters of interest are p and L = [ii, £2)- Indeed, if p = 0, then there is no 

source, while ii p > 0, the source is present at the location L. 

In the next section, we describe the Bayesian approach that will be used for model 

selection. Then the computation and algorithm will be explained. 

We use a Bayes factor approach to select between the two models Mi and M2 in 

question. The Bayes factor is defined as the ratio of the prior and posterior odds: 

pr(Mi)/pr(M2) 



BF 



pr(Mi|r)/pr(M2|F) 
pr(Mi)pr(M2)pr(F|M2) pr(F|M2) 



pr(M2)pr(Mi)pr(F|Mi) pr(r|Mi) 
This number serves as an indicator of which of the models Mi and M2 is more supported 

by the data. If BF > 1, this indicates M2 being more strongly supported by the data. 

Otherwise, Mi is more strongly supported. Furthermore, the magnitude of the Bayes 

factor is a measure of how strong the evidence is for or against Mi. According to 



Kass and Raftery 11 , when the Bayes factor exceeds 3, 20 and 150, one can say that, 
correspondingly, a positive, strong, and overwhelming evidence exists that a source is 
present. See 10 Appendix B], [6] and ^ for further interpretation of Bayes factors. 



We thus need the marginal distributions pr(y|Mj) to be calculated for each 
candidate model Mj,j = 1,2. 

Under the null model Mi, in which there is no source, one concludes that the 
corresponding marginal probability density of Y is: 

pr(f |Mi) = (47r)-'^. 

When there is a source, we denote by \1' the sample space of parameters under M2. 



The points ijj from this space is the triples ip = {p, r, u) described in Section 2.2.1 Then 
the marginal probability of Y under M2 is 

pr(F|M2)= [pT{Y,tlj\M2)f{4j\M2)d4;, 



a quantity that cannot be computed exphcitly. 
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The usual Monte Carlo method of computation is as follows. Suppose we have 
k = 1,---,K samples ip^'''^ = [p^^\r^^) ^u'^^)') from the posterior distribution. The 
marginal distribution pr(y|M2) can then be estimated as 



K 



pr(r|M2) ^ {ir-i^pr(r|^W,M2)-^}-\ 



k=l 



i.e., the harmonic mean of the likelihoods pr(y|^,M2) (see, e.g. 11 ). Given this, 
the Bayes factor is calculated as 

BF = pr(F|M2)/pr(F|Mi) (6) 

= "^"'Ef=i {(1 -p('))"-''(')(p(')rf-' + 1 -pW)^W}" 

where J{k) = J2^=i^ (Ir*^^) cos(m*^'''^ — 6i) — Si\ < d). See 11 and 15 for more details 
about calculation of Bayes factors. 

Computational details of the MCMC implementations and calculation of Bayes 
factors are shown in the following subsections. 

3.1. Computation 



With the model and prior in Section 2^, the posterior distribution is not straightforward 
to sample from. Thus, the Markov Chain Monte Carlo method is used to simulate the 
parameters from the posterior. 

Standard implementation of the Gibbs sampler in this problem will not work, 
since we discover, as Figure. [2] shows, that the posterior distributions are extremely 
multimodal. The reason for this multimodality is clear. Indeed, what the algorithm 
essentially does is to look at concentrations of trajectories at different locations. If the 
threshold is set too low, as in the left part of Figure. [21 one expects to find (and indeed 
finds) such concentrations in quite a few places. 

To overcome this problem, after some experimentation we adopted a parallel 
tempering method in order to improve mixing of simulations from this multimodal 
distribution. 

We describe now the algorithm to sample from f{ilj\Y,M2) in more details. The 



Detecting a small emission source 



9 



Figure 2. Snapshots of f{£i,£2\p,Y)- The data set contains 200,000 samples. The 
true emission rate is p = 0.001, and the source is located at (0.3, 0.6). The left figure 
is conditioned at p ~ 0.0002; the right one assumes p = 0.001. The multimodality 
illustrates the difficulty of MCMC sampling in this problem. 





reader interested only in the results of the implementation, can skip the following sub- 
sections and move directly to Section |4} 



3.1.1. Implementing the Parallel Tempering Algorithm One can find discussions of 



tempering algorithms in l8j and 13 



We run A^ parallel chains, each with equilibrium fi{x) oc f{xY^'^\ where /(x) is the 
target posterior distribution f{p,r,u\Y,M2) and Tj is a given temperature level. The 
temperature ladder Ti > ■ ■ ■ > Tjv = 1 plays the most important role in the algorithm, 
and is constructed in the following (trial-and-error) manner. We decide first the highest 
temperature such that a single MCMC run (e.g. using Metropolis-Hastings algorithm) 
at that temperature can explore the whole sample space easily (e.g. the acceptance 
rate of MH is about 90%). Then the next temperature level is chosen such that the 
rate of exchanging samples with the chain at previous temperature is moderate (e.g. 
20%). We have found in numerical experiments that A^ = 6 works well, with highest 
temperature being 5 and exponentially decreasing to 1. In our MCMC simulation, the 
Gibbs algorithm is implemented at each chain, and all chains start with random values. 






Let a;*-*^ = (x)"', ■ ■ ■ , x^/) denote the current population of samples from N chains. 



During each iteration of MCMC, the following steps of mutation and exchange are 
implemented. 
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Mutation : Update x\ to xl by the Gibbs sampler, for z = 1, ■ ■ ■ , A^. Details are 
shown in the subsection I3.1.2I 
Exchange: Starting from the first chain, try to swap with neighbors as follows, 

• For i = 1,N, exchange with the only neighbor; for i = 2,- ■ ■ ,N — 1, exchange with 
the two neighbors with equal probability. 

• Then accept the exchange of states i and j with probability 

min{l,exp([log/(x(*+^))-log/(a;f+^))] . [Tr' - Tf^])} ■ 

The chain with T/v = 1, which has the target posterior distribution as equilibrium, is 
used in the harmonic mean estimate of Bayes factors. 

3.1.2. Implementing the Gibbs Sampler In each Gibbs update, the target distribution 
is one of the /j's. In the following context, a function h{x\-) refers to the full conditional 
distribution of X, given all the other unknown variables. Notice that the unknown 
additive constants in the logarithm of a distribution do not affect Gibbs sampler. 
Joint Distribution in Each Iteration : It is easily seen that the joint posterior distribution 
of {p, r, u) is given as 

\og{Mp,r,u\Y)} = TrHog{f{p,r,u\Y)} + C 

= T-'[log{r) + \og{f{Y\p,r,u)}] + C. 

Updating p: It is easily seen that 

log{/i(p|-)} = T-'{{n - J)log(l -p) + J\og{pd-' + l-p)} + Cp. 
To update p for each fi{p\-), the following steps are taken. 

• Compute the full conditional distribution fi{pj\-) = tcj, ioi j = 1,- ■ ■ ,h. 

• Form Uj = T!'j/J2k=i'^k, for j = 1, ■ ■ ■ , h. 

• Draw from the vector (pi, ■ ■ ■ ,ph) with probabilities (cui, ■ ■ ■ , Uh)- 

Updating the Radius Component r in the Polar Coordinates: It is easily seen that 
log{/i(r|.)} = Tr'{J\og{pd-' + I - p) - J\og{l - p) + log(r)} + C,- 
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To update r, the Metropolis-Hastings algorithm is implemented. The proposal is a 
normal distribution N(rcurr, CTprop.r) truncated on the interval [0, 1], where the mean rcurr 
is the current value and the standard deviation o"prop,r is a constant. 
Updating the Angle Component u in the Polar Coordinates: It is easily seen that 

log{fi{u\-)} = Tr\j{\og{pd-^ + 1-V)- log(l -V)} + Cu. 

To update m, again the Metropolis-Hastings algorithm is implemented. Notice that the 
target function is periodic in -u, and u is restricted to [0,27r]. To ease movement across 
boundary, we propose a value Wprop from the normal distribution N(Mcurr,o"prop,u)) where 
the mean Wcurr is the current value and standard deviation o"prop,u is a constant. If 
^prop > 27r or Mprop < 0, then the candidate is reset to be Mprop,actuai = Mpropmod(27r). 

3.2. Algorithm Summary 

Our Bayesian approach could be summarized as follows. Given a particular dataset Y , 
denote the posterior by f{il)\Y.,M2) as in equation ([5|, 

(i) decide the hyperparameters Op and hp (bp > Op > 0) based on prior knowledge of p, 
and the grid ap = pi < ■ ■ ■ < ph = bp according to desired precision in p, 

(ii) decide the temperature ladder Ti > ■ ■ • > Tjv = 1 and define /«('?/') = 
{/(V^|y,M2)}i/^%for^ = l,---,iV, 

(iii) assign random initial values ipi,o = {Pio,f^io,Uio) to each MCMC chain; set t = 0, 



(iv) update ipi,t to ipi^t+i by Gibbs sampler as explained in subsection 3.1.2 , and exchange 



ipi,t+i with its neighbor(s) as explained in subsection 3.1.1, ioi i = 1,- ■ ■ ,N; set 
t = t + l, 

(v) repeat the last step until t = K; check the convergence and mixing of {4'N,t}i^o] 
adjust the temperature ladder and repeat the above steps until the MCMC chain 
converges and mixes well, 

(vi) discard the first 20% of the sequence {ipN,t}i^o (called burn-in), take every 10 
samples from the rest of the chain (called thinning), denote the new sequence by 
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{'ipj}, which are samples from the posterior distribution, 

(vii) calculate the Bayes factor estimator BF as in equation (6|, and the posterior sample 
mean ^ = {p, r, u) , 

(viii) conclude the presence of the source if BF > 3; the source strength (i.e. the emission 
rate) is estimated by p; if BF > 3, the location of the detected source is estimated 
by {f cos(M),r sin(M)}. 

Furthermore, the uncertainty in estimation could be summarized by other statistics such 
as sample standard deviation. 

4. Simulation Study 

We considered the situation where the size of the possible source is approximately known 
and is small compared to the size of the whole object. After choosing appropriate units, 
we assume that the object is the unit disk. The practically reasonable assumption is that 
the source radius is around 1% of the object radius, i.e. d = 0.01 (e.g., the object has 
dimension of several meters, while the source is of diameter of a few centimeters). The 
simulation is designed to examine the performance of the method at various emission 
rate levels, which are chosen as p = 0.01, p = 0.005, p = 0.001 and the case that no 
source exists, p = 0.00. We experimented with the number of detected particles being 
n = 2x 10^ and n = 5x 10^. With a fixed d, as the true emission rate p (and thus signal- 
to-noise ratio) decreases, a larger total number of all detected particles is required for 
detection of a source. This can be explained by a simple application of the Central Limit 
Theorem CLT (see, e.g. p]). The prior values [0 < pi, . . . ,ph] are assumed to be located 
near the true value of p, which in many applications is known with some uncertainty. 
At each level p, 10 simulated data sets where generated and analyzed, including also 
the case p = 0. The results with two different sample sizes are summarized in Table [T] 
and Table [2| respectively. Along with the Bayes factors, we also report the posterior 
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Table 1. Summary of Bayes factors for simulation in Section HI Sample size 
n — 2 X 10^. There are 10 simulations performed at each combination of level p 
and location, and 20 simulations at p = 0. The values reported in the table are 
minimum, median, maximum of 10 Bayes factors, and the proportion of Bayes factors 
being greater than 3. In the last column is the median of pr(p ~ 0\Y) calculated from 
the 10 data sets. 
p Location Min Med Max Prop>3 pr(p = 0\Y) 

0.01 (0.6,0.3) 3.7x10282 

0.01 (0.96,-0.1) 6.7x10297 

0.005 (0.6,0.3) 2.8x10^5 

0.005 (0.96,-0.1) 3.5x10^6 

0.001 (0.6, 0.3) 0.3299 

0.001 (0.96, -0.1) 0.1032 

n/a 0.36 



probability that there is no sources, namely 
pr(p = 0|F) =pr(Mi|F) 



Inf 


Inf 


1 





Inf 


Inf 


1 





3.1x10*^9 


2.4x10^9 


1 


3.1x10-9" 


1.2x10*6 


4.6xl0i"s 


1 


8.1x10-87 


11.2 


7.7x10^ 


0.6 


8.1x10-2 


1.4 


2.1x10^ 


0.4 


0.4051 


1.06 


24.73 


0.1 


0.48 



pr(Mi)pr(y|M: 



-^ = (l + BF)-\ 

pr(Mi)pr(y |Mi) + pr(M2)pr(r jMs) 

where BF refers to the Bayes factor. 

One can note from the results that if p = 0.005 or p = 0.01, much smaller sample 
sizes are sufficient to detect the source. As the level decreases, say for p = 0.0005, much 
larger sample sizes are required. In particular, the rows with p = 0.001 of the Table [T] 
show that sensitivity is not too high. The reason is that the number of detected particles, 
n = 2 X 10^, is not high enough. The next table shows significant improvements when 
the number of particles is increased. 

The results in Table |2] clearly show very high sensitivity to the presence of a source, 
which is indicated by the overall large values of BF. Furthermore, the location, if the 
source is present, can be also found with high accuracy, as shown in Figure [3) 
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Table 2. Repeat of the Tabic [l] with 5 x 10^ samples. 

p Location Min Med Max Prop>3 pr(p = 0|y) 



0.01 


(0.6,0.3) 


Inf 


bif 


bif 









0.01 


(0.96, -0.1) 


Inf 


Inf 


Inf 









0.005 


(0.6,0.3) 


1.02x10^19 


9.81x10^35 


7.77x10275 




6.48x10" 


-236 


0.005 


(0.96, -0.1) 


2.00x10^" 


4.13x10^23 


1.91x10^55 




3.35x10- 


-223 


0.001 


(0.6,0.3) 


523.76 


1.41x10^ 


5.44x10^^ 




7.32x10- 


-10 


0.001 


(0.96, -0.1) 


41.68 


2.35x10^ 


1.60x10^2 




4.96x10- 


-10 





n/a 


0.65 


1.16 


2.23 





0.46 





5. Concluding Remarks 

• The results of this study show that Bayesian methods can be successfully used for 
detection of low emission small sources in the cases of realistic parameters. High 
sensitivity can be achieved if the observation time (and thus total count of particles 
detected) is sufficiently large. 

• The assumption that the detector can determine the directional information is 
crucial. We believe that otherwise detection with such low values of SNR (signal- 
to-noise ratio) would be impossible. The determination of the incoming direction is 
usually achieved by detector collimation. This is not an option with the extremely 
low SNR levels, since it would most probably eliminate completely the useful signal. 
However, there exist the so called Compton type cameras for detecting 7-photons 
(e.g., [1,12 and references therein), as well as their analogs (although based upon a 



somewhat different physics) for neutron detection 14,17 . These cameras do not use 
collimation, but can determine some less precise directional information. Namely, 
the camera is able to provide a (hollow) cone of possible directions of the incoming 
particle. Although this is a less precise (and highly over-determined) information, 
it is known (e.g., see fl] and references therein) how to convert the Compton type 
data into the precise directional information. Thus, the algorithms described can 
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Figure 3. The estimated location of tlie source witli n ~ 5 x 10^ sample comits is 
shown for various emitting levels, with 95% highest posterior density region. The top 
left plot is with p = 0.01; the top right plot is with p = 0.005; the bottom left plot is 
with p — 0.001; the bottom right plot is with p = 0. In the plots where a source exists, 
the location of its center is indicated by the intersection of gray dashed lines. Each 
of the above figures is plotted with the posterior sample having median Bayes factors 
among the 10 simulated cases. 



» 



k:7 



0.590 0.595 0.600 0.605 



0.590 0.595 0.600 0.605 




0.590 0.595 



0.605 0.610 




be used in conjunction with Compton type detectors. 

• It is planned to address in the future study the effects of scattering, inhomogeneous 
random background noise, Compton type cameras, as well as the 3D situation. It 
is also planned to compare the results and the computational cost of the Bayesian 
approach with the more analytic techniques of [I]. 
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